#!/usr/bin/env python
import sys
import os
import hmmer3
#import MySQLdb
import csv
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from bio.parsers import fasta as fasta_parser
#import csv
#from optparse import OptionParser
'''
Use by pdzdb_build
Locating SDR from a PDZ sequence
return PDZdb tables
#----------------------
SDR table store in sdr3.tbl file
no longer need database to read sdrs
'''
def load_sdrfile():
    holder = []
    with open("sdr3.tbl") as sdrfile:
        reader = csv.reader(sdrfile, skipinitialspace=True, delimiter = '\t')
        try:
            for row in reader:
                row = [row[0], int(row[1]), int(row[2])]
                holder.append(row)
        except csv.Error, e:
            sys.exit('file %s, line %d: %s' % (filename, reader.line_num, e))
    return holder
    
def hmmscan_call(acc, hmmdb, seq, eval):

    tmpPath = "scanHmmer.tmp"
    tmpFile = open( tmpPath, "w" )
    record = SeqRecord(Seq(seq), acc, "", "")
    SeqIO.write(record, tmpFile, "fasta")
    tmpFile.close()
    
    try:
        cmdArray = ('hmmsearch -E'+str(eval)+' '+hmmdb+' '+tmpPath+
                    '> tmp.hmmscan')
        #print cmdArray
        os.system(str(cmdArray))
    except OSError:
        os.unlink( tmpPath )
        raise OSError( "Unable to find hmmsearch" )
    hmmscan_out = open('tmp.hmmscan', 'r')
    _return = hmmer3.parse(hmmscan_out)

    try:
        os.unlink( tmpPath )        
    except OSError:
        pass
    #print _return
    return _return


def sdr_search(hmmscan, motifs_set):
    total_sdrs = []
    count_domain = 0
    for record in hmmscan:
        for hit in record:
            for domain in hit:
                #if more than 1 domain found, then return a negative flag
                #means there are more than one PDZ domain in the sequence
                #and we could not tell which binds to the peptide
                if count_domain >= 1:
                    return False
                    pass
                count_domain += 1
                
                glGf_off = domain.hmm_sequence.find('glGf') + 3
                GD_off = domain.hmm_sequence.find('GD') + 1
                
                #if there's no glGf/GD motif in the hmmalignment
                #domain.target_sequence.find will return -1
                
                if glGf_off == 2 or GD_off == 0: #motif check
                    if glGf_off == 2:
                        print "GLGF motif not found in ", hit.name
                        #raw_input()
                    if GD_off ==0:
                        print "GD motif not found in ", hit.name
                        #raw_input()
                    #return an empty list, program will jump to next sequence
                    
                else:
                    #curs.execute("""select motif, offset, pep_pos from sdr;""")
                    #mySQL here is only to get info for SDR rules
                    #need to change to a file!
                    #(motifs offset and for which position)
                    #SDRs to be identified
                    #motifs_set=curs.fetchall()
                    
                    len_ts = len(domain.target_sequence)
                    sdrs = dict() #store of sdrs for each protein/one instant?
                    
                    hmm_seq = domain.hmm_sequence
                    target_seq = domain.target_sequence
                    for motif in motifs_set:
                        #=================================================#
                        #=================================================#
                        sdr = ''
                        
                        if motif[0] == "GD": #SDR offset checks
                            try:
                                count = abs(motif[1])
                                pointer = GD_off
                                
                                while count != 0:
                                    #if there are gaps in between
                                    #do not count them
                                    
                                    if motif[1] >=0:
                                        pointer += 1
                                    else:
                                        pointer -= 1
                                    if hmm_seq[pointer] not in '.-':
                                        count -= 1
                                sdr = target_seq[pointer]
                                
                            except (IndexError):
                                #print ("ERR: SDR not in hmmscan alignment",
                                #       motif[0], motif[1])
                                #raise
                                pass

                        elif motif[0] == "GLGF": 
                            try:
                                count = abs(motif[1])
                                pointer = glGf_off
                                
                                while count != 0:
                                    #if there are gaps in between
                                    #do not count them
                                    
                                    if motif[1] >= 0:
                                        pointer += 1
                                    else:
                                        pointer -=1
                                    if hmm_seq[pointer] not in '.-':
                                        count -= 1
                                sdr = target_seq[pointer]
                                
                            except (IndexError):
                                #print ("ERR: SDR not in hmmscan alignment",
                                #       motif[0], motif[1])
                                #raise
                                pass
                        else:
                            print "ERR: Unknown motif"
                        
                        if sdr not in '.-':
                            sdrs[(motif[0],motif[1],motif[2])] = sdr
                        else:
                            pass
                    #Even not all SDRs found, still record them
                    '''
                    allsdrs = True
                    for k, v in sdrs.iteritems():
                        if v == '':
                            #allsdrs = False;
                            del sdrs[k]
                    if allsdrs:
                        total_sdrs.append(sdrs)
                    '''
                    for k in sdrs.keys():
                        if sdrs[k] == '' or sdrs[k] in '.-':
                            del sdrs[k]
                    total_sdrs.append(sdrs)
    
    return total_sdrs



def build(proteinA, proteinB, pstart, pend, seq, pep, tbl_sdr, tbl_sdrdom,
          tbl_pdzdom, tbl_peptide, tbl_pepdom):
    #dbtbl.build(protein A ID(pdz), protein B ID, pdz start, pdz end,
    #            pdz sequence, peptide, sequence, five tables)
    #hmmdb = options.hmmdb
    #mydb = MySQLdb.connect(user="root", passwd="198544", db = "pdzdb2")
    #curs = mydb.cursor()
    sdrtbl = load_sdrfile()
    
    hmmscan = hmmscan_call(proteinA, 'strucaln.hmm', seq, 0.01)

    total_sdrs = sdr_search(hmmscan, sdrtbl)
    

    if total_sdrs:
        #PDZ domain (domain_ID, acc, from, to, sequence)
        #?, proteinA, pstart, pend, seq
        domainNum = 1
        domainID = proteinA + '-' + str(domainNum)

        while (domainID in tbl_pdzdom) and (
            tbl_pdzdom[domainID] != (proteinA, pstart, pend, seq)):
            #new domainID only if proteinacc etc. are different.
            domainNum += 1
            domainID = proteinA + '-' + str(domainNum)
        else:
            if domainID not in tbl_pdzdom:
                tbl_pdzdom[domainID] = (proteinA, pstart, pend, seq)
            # if domainID in tbl => this entry already in table
            # so do nothing
            
        #SDR (motif, offset, pep_pos)
        #SDR_dom(domain_ID, motif, offset, pep_pos, residue)
        if len(total_sdrs) != 0:
            for k, v in total_sdrs[0].iteritems(): 
                tbl_sdr.add(k)
                tbl_sdrdom.add((domainID, ) + k + (v, ))
        #Peptide(seq, protein_acc)
        if proteinB == pep: #proteinB is a peptide not a protein with acc
            tbl_peptide.add((pep, 'NULL'))
        else:
            tbl_peptide.add((pep, proteinB))
        
        tbl_pepdom.add((pep, domainID))

        #PDZ (acc, taxon_id)
        #PDZ Domain (domain_ID, acc, from, to, sequence)
        return tbl_sdr, tbl_sdrdom, tbl_pdzdom, tbl_peptide, tbl_pepdom
        
    else:
        #print 'No SDR found in this'
        return False
        

if __name__ == '__main__':
    sdrfile = load_sdrfile()
    print sdrfile
#    main()
